

FINAL REPORT 
NASA GRANT NAG8-704 


//- /-3 'C/Z^ 

3/ 7 


FLUID DYNAMICS AND LOW GRAVITY EFFECTS OF 
CHEMICAL VAPOR DEPOSITION 


Period of Performance 
2/22/88 through 8/21/90 


Principal Investigators 
THOMAS A. NYCE 
FRANZ ROSENBERGER 


Center for Microgravity and Materials Research 
University of Alabama in Huntsville 
Huntsville, Alabama 35899 


Table of Contents 

1. Introduction * 

2 . Numerical Modelling * 

2.1 Two-Dimensional Modelling of Horizontal CVD Reactors and its Limitations 2 

2.2 Three-Dimensional Modelling of Horizontal Chemical Vapor Deposition 2 

2.2.1 Operation at atmospheric pressure and normal gravity 2 

2.2.2 Operation at reduced pressure or low gravity 3 

2.3 Mixed Convection Horizontal Channel 4 

3. Experimental Studies of Velocity Distributions in Mixed Convection in 

Horizontal Channel 4 

3.1 Asymmetric, unsteady flow (Re = 18.5) 5 

3.2 Symmetric, steady flow with spatial oscillations(Re = 36) 5 

3.3 Symmetric, steady flow (Re = 54) 6 

4. Summary 6 

5. References ' 

6. Presentations and Publications of Research under this Grant 7 

7. Figure Captions and Figures 8 

Attachments: Reprints of Research Publications under this Grant. 


1 


1 . Introduction 

In this Final Report we summarize the work performed during the funding period of NASA 
Grant NAG8-704, i.e. from February 22, 1988 through August 21, 1990. The objective of this 
research was to provide quantitative insight into the transport conditions and resulting growth rate 
distributions in horizontal chemical vapor deposition (CVD) reactors through a complimentary 
theoretical and experimental program. In order to quantitatively test the fidelity of our numerical 
modelling, we have simulated two systems for which detailed experimental results were available: 

(a) Metal organic chemical vapor deposition (MOCVD) of gallium arsenide from 
trimethylgallium and arsine for which detailed growth rate studies had been earned out at 
atmospheric pressure by Van de Ven et al. [1]. Encouraged by the good agreement between 
numerical and experimental results, we have also simulated the operation of a few of the above 

MOCVD cases at reduced pressure or low gravity. 

(b) Mixed convection in a horizontal channel of aspect ratio 2 which we have studied by 

laser Doppler anemometry (LDA) under this Grant. 

For the MOCVD system at atmospheric pressure we will present here only a concise 
overview of the results obtained, and append our detailed publications on this topic [2,3] to this 
report. The results of the simulations at reduced pressure or low gravity will be illustrated in some 
detail; a comprehensive treatment to be published in the Journal of Crystal Growth is in 
preparation; for title see Section 7. Similarly, the results of the mixed convection studies will be 
highlighted here and presented in full detail in a later publication in the International Journal of 
Heat and Mass Transfer; see Section 7. 

2. Numerical Modelling 

Figure 1 presents the basic reactor and channel geometries that we have modelled 
numerically. We assume that a fully developed flow of, respectively, either a reactive gas mixture 
in a carrier gas (MOCVD, Fig. la) or a monocomponent gas (mixed convection, Fig. lb) flows 
from an isothermal section through a bottom heated part of same cross-section and leaves the space 
of interest through an isothermal exit section. Whereas in the MOCVD cases we have considered 
both horizontal and tilted susceptors, i.e. 0 = 0 and 0 > 0, all mixed convection simulations were 
for 0 = 0 . However, in order to explore the importance of non-perfect channel alignment for 
symmetry breaking that we observed in the experimental part of the mixed convection studies, we 
have also simulated a few cases with small tilts of the horizontal channel about its x-axis. 

For details of the assumptions (boundary conditions, kinetics, etc.) and numerical schemes 
used in the 2D and 3D modelling the reader is refereed to the attached [2,3]. The most important 
results obtained in the simulations are as follows. 
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2.1 Two-Dimensional Modelling of Horizontal CVD Reactors and its Limitations [2] 

Three models of increasing complexity were studied. It is shown that 2D models can 

produce realistic predictions only for reactors with large width-to-heigth (aspect) ratios, that are 
operated at subcritical Rayleigh number, i.e. under dominant forced flow conditions. But even 
then, depending on the carrier gas, thermal (Soret) diffusion must be included in the model. In 
comparison to pure Fickian diffusion, thermal diffusion decreases the growth at the beginning of 
the susceptor and decreases the growth rate in the downstream region. Furthermore, velocity 
corrections for finite aspect ratios must be made, and buoyancy effects can be significant in the 
entrance region. 

Under the above provisions and for hydrogen as earner gas, good agreement with 
experimental results [1] was obtained for the 2D growth rate distributions (i.e. in the axial 
symmetry plane) of GaAs throughout the susceptor. With nitrogen as carrier gas, poor agreement 
between model results and experimental findings resulted for the beginning of the susceptor , 
whereas reasonably good agreement was found for downstream regions. This was interpreted in 
terms of upstream turbulence which can be associated with the relatively high Reynolds numbers 
of the nitrogen cases, and which emphasizes the importance of proper design of the entrance flow 
section for layer uniformity. 

Our simulations show also that the results for growth distributions are relatively insensitive to 
axial diffusion, temperature dependence of transport properties (versus averaged properties at an 
average temperature), and an assumed Poiseuille velocity profile (versus a rigorous solution to the 
velocity field). 

Yet, boundary layer approximations are prone to give unrealistic results at low the Reynolds 
numbers typical for CVD operations. In contrast to the above cross-over in growth rates obtained 
from models with and without thermal diffusion, the inclusion of thermal diffusion into boundary 
layer models can only lead to either a reduction or increase in growth rate across the whole 
susceptor. 

2.2 Three-Dimensional Modelling of Horizontal Chemical Vapor Deposition 
2.2.1 Operation at atmospheric pressure and normal gravity [3] 

Full 3D, 3D parabolic and 2D solutions for the steady state transport and the resulting growth 
rate distributions were obtained for the reactors of small and large cross-sectional aspect ratios used 
in the experimental studies of [1]. Furthermore, the effects of tilting the susceptor were 
investigated for various input flow rates. We found that with light carrier gases, thermal (Soret) 
diffusion leads to more uniform growth rates in axial and cross-wise direction. Furthermore, 
depending on the aspect ratio and thermal boundary conditions on the sidewalls, buoyancy-driven 
3D flow effects can greatly influence the growth rate distribution throughout the reactor, under 
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conditions judged stable against natural convection rolls when using only criteria based on vertical 
temperature gradients (i.e. "subcritical Rayleigh numbers"). The modelling results emphasize the 
importance of the proper design of the lateral thermal boundary conditions for obtaining layers of 
uniform thickness. This study also shows that the much higher computational costs associated with 
a full 3D model, compared to 3D parabolic solutions, are justifiable only on small aspect ratio 
reactors operated at low flow rates. 

2.2.2 Operation at reduced pressure or low gravity 

Van de Ven et al's experiments [1] were conducted at atmospheric pressure and normal 
gravity, which we will refer to in the following as "standard conditions". In order to illustrate the 
consequences for layer thickness uniformity of reducing the operation pressure or the gravitational 
field, we have also simulated select 3-D cases of [3] for correspondingly changed conditions. 
Specifically, the following MOCVD cases were treated (for convenience we use here the same case 
numbers as the ones used for standard conditions, see Table 2 in [3], reprint attached): 

(1) Case 2 (carrier gas hydrogen, aspect ratio 2.8, Reynolds number 2.6, total pressure 1 atm) at 
10 -5 of normal gravity. 

Figure 2 presents the growth rate distribution in the vertical midplane for this low gravity case. 
Comparison with Fig. 8 in [3] reveals that in the midplane the layer uniformity is improved only 
slightly. Across the susceptor (i.e. normal to the main flow), however, the reduction of gravity 
leads to a drastic improvement in layer uniformity, as can be seen from the normalized growth rate 
profiles depicted in Fig. 3; for convenience we are reproducing the corresponding 1-g case (Fig. 
1 1 of [3]) in Fig. 4. 

(2) Case 6 (carrier gas nitrogen, aspect ratio 6.3, Reynolds number 5, normal gravity, partial 
pressures of reactants unchanged) at a total pressure of 0.1 atm. 

Figure 5 shows the calculated growth rate distributions for 0.1 and 1 atm, the latter 
corresponding to the heavy curve in Fig. 25 of [3]. By lowering the total pressure, the axial 
growth rate uniformity has become worse, compared to the case where the pressure is at latm. The 
plots of the normalized transverse growth rate distributions (i.e. across the susceptor) at different 
axial positions in Fig. 6, when compared to the corresponding 1 atm case depicted in Fig. 26 of 

(3) , show that for this particular case there is no advantage in operating at reduced pressure. 

(3) Case 6 (carrier gas nitrogen, aspect ratio 6.3, Reynolds number 50, total pressure 1 atm) at 

10* 3 * 5 of normal gravity. 

Figure 7, with normalized transverse growth rate distributions at various axial (downstream) 
positions, when compared to Fig. 26 of [3], shows that operation at low gravity of this specific 
case brings no advantage. 
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(4) Case 7 (carrier gas hydrogen, tilted susceptor, aspect ratio at beginning of susceptor 6.3, 
total pressure 1 atm) at 10' 5 of normal gravity. 

Comparison of Figs. 8 and 9 shows that the transverse layer uniformity of this tilted susceptor 
case can be significantly improved by operation at low gravity. 

2.3 Mixed Convection in Horizontal Channel 

In synergistic work with the LDA studies presented in Section 3, we have numerically 
modelled mixed convection in a horizontal channel of aspect ratio 2 for a variety of conditions 
closely resembling the experimental runs. The specific combinations Reynolds numbers, thermal 
boundary conditions on the sidewalls, tilt angle about the x-axis (see Fig. lb) and numerical 
scheme used are listed in Table 1. The system of conservation equations was the same as in [3], 
except for neglect of the diffusion-advection (species conservation) equation due to the 
monocomponent nature of the flow considered, and of the use of the Boussinesq approximation 
due to the small temperature difference between the hot and cold plates The meshs used had, 
respectively, 81, 30 and 15 grid points in x-, y -and z-directions, for the Re=18.75, Re=36 and 
Re=54 cases. The Rayleigh number in all simulations was 22,200, corresponding to the actual 
experimental conditions. The most important simulation results will be shown together with the 
experimental data. 

3. Experimental Studies of Velocity Distributions in Mixed Convection in a 

Horizontal Channel 

The experimental arrangement and conditions in these studies were identical to those 
employed in [4,5], except for the inner height and width of the channel, that were 2.5 and 5 cm, 
respectively. The classical hydrodynamic entrance length for a rectangular channel of this size is 
12 cm (for our maximum Re of 54) and the thermal entrance length is 8.5 cm. An isothermal 
entrance section of 40 cm in length insures that a hydrodynamically fully-developed flow was 
attained and that the temperature was uniform (at T c ) before the differentially heated section was 
reached. This was experimentally verified as well. The heated section was 81 centimeters in 
length, much longer than the thermal entrance length given (for Re=54) in the literature for 
hydrodynamically developed flows in differentially heated channels with high aspect ratios [6]. In 
all runs we used a 15°C temperature difference between the upper, hotter plate and the lower, 
colder plate of the channel. The average temperature, at which the transport coefficients of the gas 
(nitrogen) were evaluated, was 27.5°C. This resulted in a Rayleigh number, based on the channel 
height, of 22,200. Since the channel's side walls consisted of Plexiglas, the experimental thermal 
boundary conditions, relative to the thermal properties of the gas, were between adiabatic and 
conducting. This will be important for the comparison between experimental and numerical results. 
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In the following we will briefly summarize the experimental results according to the main flow 
characteristics observed. 

3.1 Asymmetric, unsteady flow (Re = 18.5) 

The main finding in this lowest Reynolds number case (average flow rate 1.18 cm/s) was an 
interesting asymmetry about the mid-width plane (y = 2.5 cm), in the flow behavior. Figure 10 
shows the evolution of u-component (axial) velocity profiles with increasing axial distance x from 
the leading edge of the heated plate, measured in the half-height plane (z = 1.25 cm). Note that 
even at x = 0, the profile is already significantly deformed as compared to the upstream Poiseuille 
profile approaching the bottom heated section. This deformation is due to a buoyancy-driven 
recirculation roll that is superimposed on the forced flow (see also Fig. 9 in [3]) and that, due to 
the low Re, actually leads to a slight flow reversal even at half height between 3.5 cm < x < 5 
cm. Most importantly, however, the maxima in the axial velocity, that arise from the two 
longitudinal rolls in the AR = 2 channel (see also [4,5]) alternate in magnitude with increasing x. 
This spatial unsteadiness gives way to an unsteady flow behavior for x > 12 cm. 

A detailed comparison of experimental and numerical results will be presented in a 
forthcoming publication. However, it is interesting to note, that the runs with perfectly horizontally 
oriented channels did not show any asymmetry about y = 2.5. Yet, simulations for a channel with 
a tilt about the x-axis as small as 2° yielded such asymmetries. This insight is particularly 
important, since in the CVD practice such small misalignments will always be present, and thus 
lateral asymmetries in the growth rates will be difficult to avoid at low flow rates. 

3.2 Symmetric, steady flow with spatial oscillations (Re = 36) 

Due to the higher average inlet flow velocity of 2.3 cm/s in this medium Reynolds number 
case, the deformation of the Poiseuille profile around x = 0 is much weaker than in the case above. 
This can be seen in the evolution of the u-component velocity profiles at z = 1.25 cm presented in 
Figs. 11-15. These figures contain also numerical solutions for these profiles obtained for adiabatic 
side walls. Note that the experimental data reflect a more rapid development (i.e. at lower x-values) 
of the flow-deforming longitudinal rolls than the numerical solutions. Preliminary numerical results 
.for conducting side walls show an even more rapid evolution of these side lobes in the u(y) 
profiles. This nicely illustrates the importance of the proper thermal boundary in modelling. At 
larger down-stream distances, however, where the mixed convection flow is more developed, the 
proper boundary conditions become less important, as one can see from the good fit of 
experimental data and numerical results for adiabatic side walls at x = 20 cm in Fig. 15. 

A detailed analysis of the many data obtained for this medium Reynolds number case shows 
that the flow is symmetric about y = 2.5 cm and steady throughout the channel. Interestingly, 
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there is still a distinct spatial oscillation in axial direction in which momentum is periodically 
transferred between the two maxima and the minimum (at y = 2.5 cm) in u(y). Figure 16 shows 
that the u(y) minima are spaced with a period of 5 cm. 

3.3 Symmetric, steady flow (Re - 54) 

As reflected by the u(y) profiles of Figs. 17, the flow in this case with the highest average 
inlet flow rate of 3.5 cm/s is symmetric about the mid-width plane and shows no spatial 
oscillations in x-direction. While this is a significant difference to the medium Re case, one can not 
exclude the possibility that such oscillations may still develop in channels longer than that used in 
our experiments. A comparison of these results with numerical solutions will be presented in the 
forthcoming detailed publication 

4 . Summary 

Based on the comparison between experimental data and numerical results for the growth of 
GaAs from TMGa, we have shown that 3D simulations are necessary to simulate rectangular CVD 
reactors even when operated under subcritical (Ra) conditions. The important points found through 
this study are summarized in the three attached reprints. 

The experimental studies of mixed convection in horizontal channels have shown three 
regimes of high Ra (22,200) number flows. At R<;=!8.5, the rolls develop very quickly, 
significantly modulating the axial velocity even before it reaches the beginning of the hot plate. A 
few centimeters downstream, the velocities become asymmetric about the vertical centerplane and 
at x=12 cm, become unsteady. These asymmetries have been predicted theoretically [7], but 
experimental evidence has not been published prior to this work. At Re=36, the axial velocity is 
only slightly modified at x=0. Although the flow remains steady and symmetric about the vertical 
centerplane, there is a small spatial oscillation in the velocities over the length of the channel. The 
period of this oscillation was around 5 centimeters. At Re=54, the longitudinal rolls developed 
smoothly over a length of 30 cm, with no asymmetries, unsteadiness, or spatial oscillations. 

Comparison of numerical simulations of these flows to experiments has revealed the 
importance and difficulty of setting proper thermal boundary conditions on the sidewalls. 
Calculated flows and experimentally measured flows showed very similar profiles, but at different 
axial locations, with the rolls developing more rapidly in the experiments. This is directly 
attributable to partially conducting sidewalls of the apparatus being hotter in the entrance section 
than the adiabatic walls of the simulations. A thorough comparison of the experimental data and 
numerical results for a variety of sidewall boundary conditions is in preparation. 
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7. Figure Captions and Figures 

Fig. 1. Definition sketch for (a) 3D CVD model configuration, with isothermal entrance 
section, bottom heated reaction section (susceptor), and isothermal exit section; (b) 
3D channel for mixed convection studies with isothermal entrance section followed 
by bottom heated section. 

Fig. 2. Growth rate distribution for case 2 in the vertical midplane of the reactor predicted 
from the full 3D model at 10' 5 g of normal gravity. 

Fig. 3. Normalized lateral growth rate distributions at 10* 5 g for case 2 predicted from the 
full 3D model at various axial z locations for half of the susceptor, (1) z = 2.5 cm; (2) 
z = 4 cm; (3) z = 8 cm; (4) z = 10 cm; (5) z = 13 cm; (6) z = 16 cm. 

Fig. 4. Normalized lateral growth rate distributions at normal gravity for case 2 predicted 
from the full 3D model at various axial z locations for half of the susceptor. Axial 
locations as in fig. 3. 

Fig. 5. Growth rate distributions for case 6 in the vertical midplane of the reactor predicted 
from the full 3-D model for operation at 1 and 0.1 atm, respectively. 

Fig. 6. Normalized transverse growth rate distributions at 0.1 atm and normal gravity for 
case 6 predicted from the full 3-D model at various axial z locations for half of the 
susceptor. Axial locations as in fig. 3. 

Fig. 7. Normalized transverse growth rate distributions at 1 atm and 10' 5 g for case 6 
predicted from the full 3-D model at various axial z locations for half of the susceptor. 
Axial locations as in fig. 3. 

Fig. 8. Normalized transverse growth rate distributions from the 3D parabolic model at 1 atm 
and normal gravity for case 7 at various axial z locations for half of the susceptor. 
Axial locations as in fig. 3. 

Fig. 9. Normalized transverse growth rate distributions from the 3D parabolic model at 1 atm 
and 10‘ 5 g for case 7 at various axial z locations for half of the susceptor. Axial 
locations as in fig. 3. 
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Fig. 

10a,b. 

Fig. 

11. 

Fig. 

12. 

Fig. 

13. 

Fig. 

14. 

Fig. 

15. 

Fig. 

16. 

Fig. 

17a,b. 


Mixed convection at Re = 18.5. Axial velocity as a function of y measured at center- 
height (z = 1.25 cm) for various axial locations. 

Mixed convection at Re = 36. Comparison of experimental and numerical axial 
velocities as a function of y at center-height (z = 1.25 cm) and x = -10 cm. 

Mixed convection at Re = 36. Comparison of experimental and numerical axial 
velocities as a function of y [ i.e. u(y)] at center-height (z = 1.25 cm) and x = 0 cm. 
Mixed convection at Re = 36. Comparison of experimental u(y) at z = 1.25 cm and x 
= 2 cm with numerical velocities. 

Mixed convection at Re = 36. Comparison of experimental u(y) at z = 1.25 cm and x 
= 5 cm with numerical velocities. 

Mixed convection at Re = 36. Comparison of experimental u(y) at z = 1.25 cm and x 
= 2 cm with numerical velocities. 

Axial velocity as afunction of axial distance x, measured along centerline (z = 1 .25 
cm, y = 2.5 cm) for Re = 36. 

Mixed convection at Re = 54. Axial velocity as a function of y measured at center- 
height (z = 1.25 cm) for various axial locations. 
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FIG. 17b 




